# define which candidate interface TFs are reachable through expressed signalling paths from source nodes
# (remove not-expressed proteins from pathways, collect TFs still connected)

library(igraph)

args = commandArgs(TRUE)
cutoff = args[1]

# Metacore interactions
obj_table = read.table("../data/table-pathwayObjects-Metacore.txt",sep="\t",quote="")
TF_table = read.table("../data/tableTFpath-Metacore.txt",sep="\t",quote="")

# read 1-(booleanization p-value)
pexpr = read.table("bool_MetacoreObjects.txt",stringsAsFactors=FALSE,row.names=1,sep="\t",quote="")
# define which nodes are not expressed
not_expr = rownames(pexpr)[which(pexpr[,2] < cutoff)]

connected_TFs = list()
for (pathway in colnames(obj_table)){
	# prepare pathway network
	pathway_graph = readRDS(paste0("../data/graph_str-",sub(".","-",pathway,fixed=TRUE),".Rdata"))
	if ("" %in% V(pathway_graph)$name) pathway_graph = delete_vertices(pathway_graph,"")
	# identify source nodes for this pathway
	sources = names(which(degree(pathway_graph,mode="in",loops=FALSE)==0))
	# add a source dummy node upstream of all source nodes
	pathway_graph = add_vertices(pathway_graph,1,name="dummy")
	pathway_graph = add_edges(pathway_graph,edges = as.vector(t(cbind("dummy",sources))))
	# remove the non expressed signalling nodes
	to_remove = rownames(obj_table[not_expr,])[which(obj_table[not_expr,pathway]==1)]
	expr_pathway_graph = induced_subgraph(pathway_graph,setdiff(V(pathway_graph)$name,to_remove))
	# select the TFs
	pathway_TFs = rownames(TF_table)[which(TF_table[,pathway]==1)]
	pathway_TFs = pathway_TFs[pathway_TFs %in% V(expr_pathway_graph)$name]
	#collect TFs that are disconnected from dummy node or not expressed themselves
	missing_TFs = pathway_TFs[which(pathway_TFs %in% to_remove)]
	d = distances(expr_pathway_graph,"dummy",setdiff(pathway_TFs,missing_TFs),"out")
	missing_TFs = c(missing_TFs,setdiff(pathway_TFs,missing_TFs)[which(!is.finite(d))])
	# define which TFs stay
	remaining_TFs = setdiff(pathway_TFs,missing_TFs)
	connected_TFs[[pathway]] = remaining_TFs
}

# transform into table
connected = unique(unlist(connected_TFs))
connected_table = matrix(0,nrow = length(connected),ncol = length(connected_TFs),dimnames=list(connected,names(connected_TFs)))
for (x in seq_along(connected_TFs)){
	connected_table[connected_TFs[[x]],x] = 1
}
write.table(connected_table,sprintf("table-reachableTFs-cutoff%s.txt",cutoff),sep="\t",quote=FALSE)

